A GEANT-based study of atmospheric neutrino oscillation parameters at INO 
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We have studied the dependence of the allowed space of the atmospheric neutrino oscillation parameters on 
the time of exposure for a magnetized Iron CALorimeter (ICAL) detector at the India-based Neutrino Obser- 
vatory (INO). We have performed a Monte Carlo simulation for a 50 kTon ICAL detector generating events by 
the neutrino generator NUANCE and simulating the detector response by GEANT. A chi-square analysis for 
the ratio of the up-going and down-going neutrinos as a function of L/E is performed and the allowed regions 
at 90% and 99% CL are displayed. These results are found to be better than the current experimental results of 
MINOS and Super-K. The possibilities of further improvement have also been discussed. 

PACS numbers: 14.60.Pq, 96.40.Tv 
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I. INTRODUCTION 

The evidence of neutrino masses and their mixing has 
brought neutrino physics into centre stage of particle physics. 
The neutrino mass eigenvalues and the Pontecorvo, Maki, 
Nakagawa, Sakata (PMNS) mixing matrix [3, 4] connecting 
the mass to the flavor basis provides a natural framework for 
handling three active neutrinos. 

The present information on the neutrino mass-squared dif- 
ferences and mixing angles are the following: From atmo- 
spheric neutrino detection one gets the best-fit values with 3o 

error jAm^jl - '^■^-ol 1^"^ ^^2, sin^Gis ^ 0.5to[n while 
solar neutrinos tell us Amji ~ 7.9 x 10"^ eV^, sin^9i2 — 0.30 
Jit]. Here we define Am? = m? — m?. 

At the moment, the sign of Amjj is not known. The posi- 
tive/negative value of this quantity denotes the direct/inverted 
mass ordering. The two large mixing angles and the mass 
squared differences may permit measurement of CP-violation 
in the lepton sector, if the third mixing angle, 613, and the CP 
phase, 5, are not too small. The current bound on the former 
is sin^ 9i3 < 0.05 (3a) 10. Et while 5 is unconstrained. 

Thus the determination of mass hierarchy and the measure- 
ment of oscillation parameters with high precision are of ut- 
most importance. Also of importance is observing a full 
oscillation cycle to convincingly establish that it is truly neu- 
trino oscillation which is at play. Most experiments observe 
the depletion part but not the regeneration part of the cycle. A 
reanalysis of old Super-Kamiokande (Super-K) data claimed 
to observe this 18]. However, a reconfirmation of this with bet- 
ter statistics is much awaited. The mixing angle and the mass 
squared difference for the atmospheric sector should also be 
measured more accurately. The sensitivity of the measure- 



ment of a particular parameter depends crucially on the ranges 
of neutrino energy and path length traversed from the source 
to the detector These ranges can be set in case of neutrino 
beams from artificial sources like nuclear reactors (energy ~ 
MeV) and accelerators (energy ^ GeV). Neutrinos with en- 
ergy ^ MeV (GeV) can also be obtained from natural sources 
like the sun (the atmosphere). Unlike typical accelerator or 
reactor neutrinos, the spectrum of atmospheric neutrinos cov- 
ers many decades of energy (E ^ 100 MeV - few hundred 
GeV) with comparable interaction rate and baseline (L ~ 10 
km - 12800 km). Since the oscillation probability depends 
mainly on L/E which varies in a wide range for atmospheric 
neutrinos, the measurement of the appearance/disappearance 
probability as a function of L/E can explore its variation over 
this entire range. This advantage is partly offset, however, by 
the difficulty that the flux is less known compared to that from 
man-made sources. 

Currently around the world, there are many ongoing and 
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planned experiments: MINOS 

[12, 13J, NOvA 0113, Double Chooz Slm^UNO uJ], 
Super-K III Hyper-K SElll, OPERA El 03 etc. 
Out of these only MINOS employs a magnetic field and has a 
good charge identification capability. It is to be noted that all 
these experiments are planned in the northern hemisphere of 
the earth. 

The proposed India-based Neutrino Observatory (INO) 
ll25ll at a site close to the equator plans to use a large mag- 
netized Iron CALorimeter (ICAL) detector. The proposal is 
for an underground facility with more than 1 km overburden. 
Since the detector has a high charge identification capability 
(> 90% after selection of events as described in section V and 
70% before doing the selection) ll25ll . it has a good chance 
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of determining the neutrino mass ordering 1261 1271 l28l 12911 
and also of studying the deviation from maximahty for 623 
iSElll- 

In this work we first demonstrate through a GEANT-based 
simulation of atmospheric neutrinos that ICAL indeed is ca- 
pable of observing the full oscillation cycle. We have used 
a two flavor oscillation formalism and studied the precision 
that can be achieved for |Am^| and sin^29 at INO with atmo- 
spheric neutrinos. Though a more realistic approach would be 
the use of three flavor analysis, the smallness of the mixing 
angle 813 ensures that the two flavor approximation mimics 
the real situation reasonably well. The precision depends on 
the exposure in terms of kTon-yr, reconstruction method, and 
the selection of the events in the analysis. The paper is orga- 
nized as follows: A brief summary of the neutrino oscillation 
formalism is given in Section II. The ICAL detector at INO is 
described in Section III. In Section IV a brief account of the 
atmospheric neutrino flux that has been used in the present 
analysis has been furnished. The generation of simulated data 
at ICAL and the analysis of such data are described in Section 
V. In Section VI we present the results and precision study of 
the oscillation parameters. Finally, Section VII includes dis- 
cussions and conclusions. 



II. NEUTRINO OSCILLATION 

A neutrino flavor eigenstate |Va} (oc = e,/j,T etc.) can be 
written as a linear superposition of neutrino eigenstates |v,} 
(with definite non-degenerate masses mf) in the mass basis as 
|Va) — Lft^aiiv,) (i — 1,2,3 etc.). Here Uai are the matrix el- 
ements of the neutrino mixing matrix U. This gives rise to the 
phenomenon of neutrino flavor oscillation. The probability 
that a neutrino Vg with energy E gets converted into another 
neutrino V/ after traversing a distance L in vacuum is given by 



P(V, ^V/) = 5/,-4^Re(f/;,f/,,.t//,-t/;^>in2(1.27Am2.-) 



±2£lm(t/;,f/„t/^,f/;^.)sin(2.54Am2.-)(l) 



In the above, L is expressed in km, E in GeV and Arn^ in eV^. 
The - (+) refers to neutrinos (anti-neutrinos). 

For a two flavor scenario the above equation takes a simpli- 
fied form given by 



survival 



1 -sin2 2esin2(1.27Am2-) 



(2) 



III. THE INO DETECTOR 

The simulation has been carried out for a detector with 50 
kTon mass with dimension 48 mx 16mx 12m for ICAL 
ll25ll . The detector consists of a stack of 140 horizontal lay- 
ers of 6 cm thick iron slabs interleaved with 2.5 cm gap for 
the active detector elements. For the sake of illustration, we 
define a rectangular coordinate frame with origin at the center 
of the detector, x(3')-axis along the longest (shortest) lateral 
direction, and z-axis along the vertical direction. A magnetic 
field of strength 1 Tesla is considered along the +3'-direction. 
Resistive plate chambers (RPC) have been chosen as the ac- 
tive part of the detector The readout of the RPCs is through 
the Cu strips having 2 cm width and placed orthogonally on 
the two external sides of the detectors. This type of detector 
has good time (~ 1 ns) and spatial resolutions. 



IV. ATMOSPHERIC NEUTRINO FLUX 

The atmospheric neutrinos are produced from the interac- 
tions of the cosmic rays with earth's atmosphere. The knowl- 
edge of the primary spectrum of cosmic rays has been im- 
proved from the observations by BESS 11321 13311 and AMS 
ll34ll . However, a large region of parameter space has been 
unexplored and they are interpolated or extrapolated from the 
measured flux. The difficulties and uncertainties in the cal- 
culation of the neutrino flux depend on the neutrino energy 
II35II . The low energy flux is known quite well. The cosmic 
ray fluxes (< 10 GeV) are modulated by solar activity and 
the geomagnetic field through a rigidity (momentum/charge) 
cutoff. At higher neutrino energy (> 100 GeV), solar activ- 
ity and rigidity cutoff are irrelevant. There is an agreement 
within 5% among the calculations for neutrino energy below 
10 GeV though different groups used different hadronic inter- 
action models in their calculations. 

We use the neutrino interaction cross section model of NU- 



where 9 is the mixing angle of neutrinos. Herein and in the 
rest of the paper the symbol 9 and Arn^ refer to Q23 and Amgj. 



ANCE ll36ll incorporating a typical Honda flux calculated in a 
3 -dimensional scheme 13511. 



V. DATA GENERATION AND METHOD OF ANALYSIS 

The interactions of neutrinos with the detector material are 
simulated by the Monte Cai-lo method in NUANCE In 
order to study the ICAL detector response for each event we 
use another Monte Carlo code GEANT Q- The GEANT 
code uses the information of vertex position and momentum 
of the product particles obtained from the output of NUANCE 
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simulation. 

Event reconstruction: Our present analysis is based on 
the tracks generated by the muons that are produced by the 
charged current (CC) interactions of the neutrinos in the de- 
tector volume. The muons lose energy mainly due to ioniza- 
tion and radiative processes. This energy loss is proportional 
to the effective path-length which is the product of geomet- 
ric path-length and the density of the medium. This can be 
applied only for fully contained (FC) events. 

In this simulation study we do not consider any atmospheric 
muon background and the noise hits produced by the detector. 
However, we conservatively assumed hits > 6 will be required 
for the reconstruction and filtering of the muon events from 
the latter background. 

For a given triggered energy of a muon, the number of hits 
decreases when one goes from vertical to horizontal direc- 
tion since it traverses less number of active detector elements. 
This dependence on the direction is less for the effective path- 
length. 

In case of partially contained (PC) events the momentum 
has been determined from the curvature of the track at the 
vertex due to the magnetic field applied across the detector. 
Due to limitation of our PC event reconstruction algorithm 
we considered PC events with hits < 20. 

On the other hand, the hadrons produced create a shower of 
hits around the vertex of the event. This implies that for any 
particular CC event the longest track normally comes from 
muons and this can be utilized for the analysis. 

The up-going muon type neutrinos traverse larger path- 
length undergoing oscillation whereas the down-going ones 
with much shorter path-length have little chance to oscillate. 
So one can visualize our detector set up as far (near) detector 
for the up-going (down-going) neutrinos. Then the ratio of 
up-going and down-going neutrinos (up/down) will roughly 
mimic the survival probability. This up/down ratio as a func- 
tion of L/E minimizes the systematic uncertainties in flux as 
well as in cross sections. The length L for the up-going neu- 
trinos is the actual path-length traversed by them whereas for 
down-going neutrinos the reference path-length L is consid- 
ered to be that of associated up-going neutrinos with zenith 
angle (180° — 9zenith) so that the range of L/E remains the 
same for up-going and down-going neutrinos 113 811 . 



Selection of events and Resolutions: 

The L/E resolution has a complicated dependence on L and 
E. However, a few general remarks can be made here. Qual- 
itatively for a fixed energy, the L resolution worsens gradu- 
ally as we go from vertical to horizontal region and worsens 
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FIG. 1: The selection of events in E-L/E plane, which gives 
good optimization between statistics and L/E resolution. 



rapidly close to the horizon. Also for a fixed direction, the 
L/E resolution improves with increase in E. If one neglects 
totally the near horizon events (say, between zenith angle 70° 
and 110°) all the events below 200 km/GeV are lost. In our 
analysis we consider only the high energy events at near hori- 
zon and relax it gradually as we move away from the horizon. 
Quantitatively, this is taken care of by an E dependent cut of 
the form : 

L/E> aE'' 

broken into three segments as shown in fig. [ij 

Using the above cut the resolutions for E, L and L/E ob- 
tained with the atmospheric neutrino flux for the whole range 
of E^ and Lv are shown in fig. |2l 

A representative statistics for 5 year data is shown in table 
U Here the number of events with hits > 6 is considered for 
the analysis. The selected number of events is further reduced 
by imposing the above E dependent L//? cut for a better reso- 
lution. 



cut 


No. of survivin 


I events/efficiency 




FC 


efficiency 


FC+PC 


efficiency 


hit> 6 


4160 




5351 




(for reconstruction) 










E dependent L/E cut 


2089 


50.2% 


2808 


52.4% 


(improves L/E resolution) 











TABLE I: Sample number of events after cuts in 5-year data for 

Am^ = 2.3 X lO^^^eV^. 



Using the above cuts we find the up/down distribution for 
different time exposures of the ICAL detector Here our main 
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FIG. 2: The ICAL resolutions for E, L and L/E with atmospheric neutrinos for the whole range of £v and Lv. The subscript 'ex' and 
V refer to the reconstructed and true v values. 
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FIG. 3: The simulated up/down distribution at ICAL as a function of L/E for 5 years FC events with hn? = 2.3 x 10 ■'eV^. 



goal is to find how precisely one can measure Am^. A rep- 
resentative 'up/down' distribution with respect to L/E for 5 
year FC events is shown in fig. [3] Such simulated plots are 
referred to as the 'experimental up/down' distributions. 

X^-analysis: 



In the x^-analysis the 'theoretical up/down' distribution 
is obtained by taking 40 years of atmospheric un-oscillated 
charged current muon neutrino data. The oscillation proba- 
bility is then calculated from the L and E of the neutrinos for 
each event and the event is kept or rejected by throwing a ran- 
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dom number. We smear this over the whole range of L/E 
following the L/Zs -resolution function. Finally the up/down 
ratio is calculated for different L/E bins. In this process we 
are also minimizing the effects due to geomagnetism and the 
shape of the earth. 

Then a J^-fii is made with the 'experimental up/down' dis- 
tribution varying the atmospheric mass square difference Am^ 
and the mixing angle 9 in the 'theoretical up/down' data. 



VI. RESULTS 

The L/E plot in fig. |3] shows clearly a full oscillation cy- 
cle and is typical for the ICAL detector Thus ICAL is in a 
position to observe the oscillation pattern better than previous 
attempts, hke the one by Super-K [SI]. 

We show in fig. |4]the contours in the Am^ — sin^29 plane 
for 90% and 99% CL with 5 year FC (upper left) and 5 year 
FCh-PC (upper right) events for the input value of Am^ = 2.5 x 
10^'^eV^ and with 10 years FC events (lower) for Am^ — 2.7 x 
10--^eV2. 

It is noted that the extracted best-fit values gradually be- 
come close to the input value of Am^ with increase of exposure 
time. For example, the best-fit value is found to be 2.50, 2.40 
(2.50, 2.34) xlO-^eV^ for 5, 10 years FC (FCh-PC) events 
with the input 2.3 xlO^^eV^. For all these cases the best- 
fit values of sin^ 20 turns out to be 1 with the input value 1 . 
Since the FCh-PC sample contains more high energy events 
than the FC sample and the E and L resolution are better for 
PC events, the best-fit values obtained from the FC+PC anal- 
ysis are closer to the input value. With a change of the input 
Am^ from 2.5 to 2.7 x lO^^eV^ for 10 year FC samples, the 
best-fit changes to 2.46 to 2.68 x lO^^eV^. 

The position of the dip in this up/down distribution is in- 
dicative of the best-fit value of Am^ while the overall statistics 
determines the size of the allowed parameter space. Particu- 
larly, the statistics in the larger (smaller) L/E region from the 
dip determines the lower (upper) bound of Am^. For atmo- 
spheric neutrinos the statistics increases with increase of L/E 
thus resulting in a better lower bound in the contour of Am^. 
We have made a comparison of the present results and those 



obtained by the Super-K |39D and the MINOS |110|] exper- 







iments. They are shown in fig. |5] Here we plot contours at 
90% and 99% CL in Am^ — sin^ 29 plane as obtained from our 
study of 5 years FC data and those obtained from previous 
Super-K (1489 days data) and the recent MINOS data. For 
the case of Super-K, two different analyses, one with respect 
to zenith angle and the other with respect to L/E |^8^ are 
given. One can see clearly that ICAL results are far more pre- 



cise than those of MINOS and substantially better than those 
from Super-K. 



A. Precision 

We define precision (P) for a certain confidence level of a 
particular set of oscillation parameters as 

UL-LL 



UL + LL 



where 'UL' and 'LL' are the upper and lower limit of the con- 
tour respectively at the specified confidence level. 

FC analysis: The variation of precision of sin^ 9 and Am^ 
with different years of exposure is shown separately in fig. |6] 
at 90% and 99% CL. It is seen that the precision falls very 
slowly beyond ten years and that can be a useful observation 
for the future experiment. 

It is further observed that the precision gradually becomes 
worse when we increase the value of Aw^ from 2.3 x lO^^eV^. 
This is demonstrated in fig. [T] 

The reason behind this is the following. The position of 
the dip in the up/down distribution shifts towards larger val- 
ues of L/E with the decrease of the value of Arri^. The flux 
increases rapidly with decrease of energy and the statistics be- 
comes gradually high at larger L/E. 

However, we comment that there is a competition between 
statistics and L/E resolution. We see that at low value of Am^, 
say 2.1 xlO-^eV^, the precision worsens compared to that at 
2.3 x lO^^eV^. This can be improved if we choose more strin- 
gent cuts with good L/E resolution for larger values of L/E. 
So one has to optimize between the requirement of statistics 
and L/Zi-resolution, which depends mainly on the range of 
interest of Am^. 

FC+PC analysis: After the inclusion of the PC events the 
results are very similar to that obtained from FC events and 
hence are not shown separately. 



VII. DISCUSSIONS 

Simulation studies for atmospheric neutrinos at the pro- 
posed Iron Calorimeter (ICAL) detector at INO have been 
made with a goal to determine the level of precision which 
may be achieved. The oscillated atmospheric neutrino events 
for a known set of values of oscillation parameters are gener- 
ated with the event generator code (NUANCE) and the simu- 
lated signals in the detector are obtained through the detector 
simulation code (GEANT) that uses the NUANCE output as 
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FIG. 4: The contour plots in Am^ - sin^2e plane with input hn? = 2.3 x 10"-' eV^ for 5 years FC events (upper left), FC+PC event 
(upper right); the lower plot is with input Am^ = 2.7 x lO^-'eV^ for 10 years FC events. 



0.005 
0.0045 

0.004 
0.0035 

0.003 
0.0025 

0.002 
0.0015 

0.001 



ICAL (5 yr FC) CL 90% ■ 
ICAL (5 yr FC) CL 99% ■ 
MINOS CL 90% . 



SK (zenith) CL 90% 
SK (zenith) CL 99% 
SK (L/E) CL 90% 
SK (L/E) CL 99% 




0.7 



0.75 



0.8 0.85 

sin220 



0.9 



0.95 



FIG. 5: The contours at 90% and 99% CL for 5 years FC events with An? = 2.3 x lO^^eV^ at ICAL and the contours from the 
current experiments. 



its inputs. A %^ analysis of the results obtained from this sim- 
ulated GEANT output data, properly chosen using appropriate 
constraints ("cuts"), is performed for the precision studies. 
There is, however, scope for improvement of these studies. 

• The present analysis is performed only with the simu- 
lated muon signals neglecting the hadrons. The estima- 
tion of neutrino energy E and L/E is expected to im- 
prove by the inclusion of hadrons. As hadrons mainly 



produce showers instead of well defined tracks, the 
method of using tracklength is not effective to extract 
energy information (or directional information) from 
such hadron showers. A new methodology is to be de- 
veloped for this purpose and this work of incorporating 
hadrons in the analysis is in progress. 

• A full three flavor analysis can address issues like mat- 
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FIG. 6: The variation of the precision of Am^ (left) and sin^9 (right) with time of exposure for 50l<Ton ICAL with FC events for the 
input of Am^ = 2.3 X lO^^eV^. 
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FIG. 7: The variation of the precision of Am^ (left) and sin^9 (right) with the input value of Am^ for SOkTon ICAL with 10 years FC 
events. 



ter effect, mass hierarchy and the deviation from maxi- 
maUty of the atmospheric mixing. 

• For the analysis of PC events the curvature of the track 
is used for calculation of energy. In the high energy 
regime the tracks have no or negligible curvature inside 
the detector volume and hence such PC events could not 
be considered in the present analysis. 

Moreover, we find the resolutions are energy dependent and 
significantly different for neutrinos and anti-neutrinos. One 
therefore expects to obtain more precise best fit values with 



CL contours in parameter space further shrunk, if one uses 
multiple resolution functions, instead of one as used here. 
In doing so, the whole L — E plane is divided into multiple 
small segments (mesh) and separate resolution functions are 
obtained for each such segment of the mesh which is then used 
for the purpose of analysis. 
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